Size Spectrum Sensitivity
Adam A. Kemberling
9/15/2020
Log-length and log weight bins can be created using case_when()
# Add some log bins
dataBin <- dataBin %>%
mutate(
lweight = log(wmin / 1000),
llen = log(LngtClass),
l10wt = log10(wmin / 1000))
# length bins
dataBin <- dataBin %>%
mutate(
lenbin = case_when(
llen > 5 ~ 10,
llen > 4.5 ~ 9,
llen > 4 ~ 8,
llen > 3.5 ~ 7,
llen > 3 ~ 6,
llen > 2.5 ~ 5,
llen > 2 ~ 4,
llen > 1.5 ~ 3,
llen > 1 ~ 2,
llen > 0.5 ~ 1,
TRUE ~ 0))
# weight bins
dataBin <- dataBin %>%
mutate(
wtbin = case_when(
lweight > 4 ~ 8,
lweight > 2 ~ 7,
lweight > 0 ~ 6,
lweight > -2 ~ 5,
lweight > -4 ~ 4,
lweight > -6 ~ 3,
lweight > -8 ~ 2,
lweight > -10 ~ 1,
lweight < -10 ~ 0))
# weight log10 bins
dataBin <- dataBin %>%
mutate(
wt10bin = case_when(
l10wt > 2.5 ~ 20,
l10wt > 2 ~ 19,
l10wt > 1.5 ~ 18,
l10wt > 1 ~ 17,
l10wt > 0.5 ~ 16,
l10wt > 0 ~ 15,
l10wt > -0.5 ~ 14,
l10wt > -1 ~ 13,
l10wt > -1.5 ~ 12,
l10wt > -2 ~ 11,
l10wt > -2.5 ~ 10,
l10wt > -3 ~ 9,
l10wt > -3.5 ~ 8,
l10wt > -4 ~ 7,
l10wt > -4.5 ~ 6,
l10wt > -5 ~ 5,
l10wt > -5.5 ~ 4,
l10wt > -6 ~ 3,
l10wt > -6.5 ~ 2,
l10wt > -7 ~ 1,
TRUE ~ 0)) Load results from
# 5g min size cutoff results & the stratified abundance
# SOURCE: 08_nefsc_sensitivity_prep.R
min5g <- read_csv(here("data/size_spectra_results/nefsc_5grams.csv"))
min5g_strat <- read_csv(here("data/size_spectra_results/nefsc_5grams_strat.csv"))
# Load OISST Timeseries for Survey Regions
regions <- c("Georges Bank", "Gulf of Maine", "Southern New England", "Mid-Atlantic Bight")
region_temps <- map(regions, function(x){
region_name <- str_replace_all(str_to_lower(x), " ", "_")
region_path <- str_c("~/Box/NSF OKN Demo Data/oisst/regional_timeseries/nmfs_trawl_regions/OISSTv2_anom_",
region_name, ".csv")
region_ts <- read_csv(region_path)
return(region_ts) })
# Set names, append
region_temps <- region_temps %>% setNames(regions) %>% bind_rows(.id = "region")NEFSC NMFS Trawl Size Spectra Sensitivity Exploration
Goal of this analysis is to explore the sensitivity of NEFSC size spectra slopes to adjustments in minimum/maximum bodymass cutoffs.
Additional exploration will be done to examine catch composition through:
Percent of total biomass sourced from different classes of fishes
Cumulative distribution figures, examining the % of total biomass within discrete weight bins
Detailed code for the data builds and summary functions can be found here.
Minimum Bodymass Sensitivity
One way to address differences in catchability of different sized fishes is by setting a minimum biomass threshold for the size spectra analysis. Ideally the minimum biomass threshold is set so that you retain all sizes that are reliably caught by the sampling gear. You also would like to see the size-spectrum equations to not shift dramatically by small adjustments in this cutoff.
For a starting point in the sensitivity analysis I will use all the data and calculate slopes for all the data across an evenly spaced range of cutoff values similar to a grid search.
# Range of size cutoffs
size_cutoffs <- c(seq(0, 50, by = 5))
# Calculate the slope/intercepts
size_cutoff_results <- map(size_cutoffs, function(cutoff_size) {
# Filter with cutoff size
filtered_data <- dataBin %>%
filter(wmin >= cutoff_size)
# Set grouping level
grouped_data <- filtered_data %>%
mutate(group_level = "all_data") %>%
split(.$group_level)
# get SS results
group_results <- grouped_data %>%
imap_dfr(group_mle_calc) %>%
mutate(stdErr = (abs(confMin - b) + abs(confMax - b)) / (2 * 1.96),
Year = "all",
season = "all",
area = "all",
C = (b != -1 ) * (b + 1) / ( xmax^(b + 1) - xmin^(b + 1) ) + (b == -1) * 1 / ( log(xmax) - log(xmin)))
return(group_results)
})
# Set the names for the list
size_cutoff_results <- size_cutoff_results %>%
setNames(str_c(size_cutoffs, " grams")) %>%
bind_rows(.id = "cutoff_g")
# Format the Slope results so the cutoff is numeric
size_cutoff_results <- size_cutoff_results %>%
mutate(cutoff = str_sub(cutoff_g, 1, -7),
cutoff = as.numeric(cutoff)) Slope Sensitivity
Slope seems to decline in an nearly linear way as minimum bodysize cutoff increases. Thi linear relationship makes sense because you are removing the biomass of smaller fishes, meaning the remaining larger fishers take up a larger proportion of the total biomass. This makes me wonder whether the actual slope value is important, or whether the relative changes in slope may be more important. It would seem like as long as the cutoff is consistent between groups that they could be compared against one-another for inference.
In the Edwards vignettes a cutoff of 5g is used, which is referenced to a paper by Blanchard et al., 2005.
# Plot the differences in SS slopes
size_cutoff_results %>%
ggplot(aes(x = cutoff, y = b)) +
geom_point() +
geom_smooth(formula = y ~ x, method = "lm", se = TRUE, size = 1) +
stat_poly_eq(formula = y ~ x, parse = TRUE,
label.y = "top", label.x = "right") +
labs(x = "Minimum Bodymass Cutoff (g)", y = "Size Spectrum Slope (b)")
Biomass Bins
Taking some coarse breaks in
# 4 plots for how the distribution looks with those cutoffs
coarse_breaks <- c(0, 250, 500, 750, 1000)
coarse_breaks <- setNames(coarse_breaks, paste0(coarse_breaks, "g"))
# bin summaries for overall view
break_plots <- map(coarse_breaks, function(breakpoint){
bio_plot <- dataBin %>%
filter(wmin >= breakpoint) %>%
group_by(lenbin, wtbin, wt10bin) %>%
summarise(n_group = sum(Number)) %>%
mutate(mass_break = paste0("Minimum Body Mass = ", breakpoint)) %>%
ggplot(aes(wt10bin, n_group, fill = mass_break)) +
geom_col(show.legend = F) +
scale_y_log10() +
scale_fill_gmri() +
facet_wrap(~mass_break) +
labs(x = "log10 Weight Bins", y = "Abundance") +
xlim(c(0, 24))
})
break_plots$`0g`
break_plots$`250g`
break_plots$`500g`
break_plots$`750g`
break_plots$`1000g`
Biomass Allocation
Independently of the slope and intercept for each group there may be value in seeing where across the spectrum of individual biomasses the majority of bodymass resides.
Not sure if cumulative density plots of stacked barplots are more helpful so I will explore both.
# cumulative distribution data prep
# Annual Totals
totals <- dataBin %>%
group_by(Year) %>%
summarise(`Total Annual Biomass` = sum(Biomass)) %>%
ungroup()
# Size Bin Data Prep
bin_totals <- dataBin %>%
group_by(Year, `Bodymass Bin (g)`) %>%
summarise(`Total Bin Biomass` = sum(Biomass)) %>%
ungroup() %>%
left_join(totals, by = "Year") %>%
mutate(`Percent of Annual Biomass` = (`Total Bin Biomass` / `Total Annual Biomass`) * 100 )
# Species Class Bins Data Prep
class_totals <- dataBin %>%
group_by(Year, spec_class) %>%
summarise(`Total Class Biomass` = sum(Biomass)) %>%
ungroup() %>%
left_join(totals, by = "Year") %>%
mutate(`Percent of Annual Biomass` = (`Total Class Biomass` / `Total Annual Biomass`) * 100 ,
`Species Category` = spec_class)
########## Cumulative Sum Plots Prep ##################
# 1. Data Prep for cumulative sums on size bins
size_cumsums <- bin_totals %>%
arrange(Year, `Bodymass Bin (g)`) %>%
group_by(Year) %>%
mutate(`Cumulative Biomass (kg)` = cumsum(`Total Bin Biomass`) / 1000,
`Proportion of Annual Biomass` = `Cumulative Biomass (kg)` / (`Total Annual Biomass` / 1000),
`Weight Bin` = as.numeric(factor(`Bodymass Bin (g)`)),
`Weight Bin` = as.factor(`Weight Bin`))
# 2. Cumulative sums on species class
class_cumsums <- class_totals %>%
arrange(Year, `Species Category`) %>%
group_by(Year) %>%
mutate(`Cumulative Biomass (kg)` = cumsum(`Total Class Biomass`) / 1000,
`Proportion of Annual Biomass` = `Cumulative Biomass (kg)` / (`Total Annual Biomass` / 1000),
`Species Class` = as.numeric(factor(`Species Category`)),
`Species class` = as.factor(`Species Class`))
# 3. Cumulative Sums on lengths
lmin_cumsums <- dataBin %>%
group_by(Year, LngtClass) %>%
summarise(`Total Bin Biomass` = sum(Biomass)) %>%
ungroup() %>%
left_join(totals, by = "Year") %>%
mutate(`Percent of Annual Biomass (g)` = (`Total Bin Biomass` / `Total Annual Biomass`) * 100 ) %>%
arrange(Year, LngtClass) %>%
group_by(Year) %>%
mutate(`Cumulative Biomass (kg)` = cumsum(`Total Bin Biomass`) / 1000,
`Proportion of Annual Biomass` = `Cumulative Biomass (kg)` / (`Total Annual Biomass` / 1000),
`Length (cm)` = as.numeric(factor(LngtClass)))Bin Size Proportions
# Size Bins merging with totals - This one could be better as cumulative distirbution plot (s)
bin_totals %>%
ggplot(aes(x = fct_rev(Year), y = `Percent of Annual Biomass`, fill = `Bodymass Bin (g)`)) +
geom_col(color = "white", size = 0.2) +
labs(x = "", y = "Proportion of Total Annual Biomass") +
coord_flip() +
scale_fill_gmri() +
scale_y_continuous(position = "right") +
guides("fill" = guide_legend(title.position = "top", title.hjust = 0.5)) +
theme(legend.position = "bottom")
Species Class Proportions
# Class Proportions
class_totals %>%
ggplot(aes(x = fct_rev(Year), y = `Percent of Annual Biomass`, fill = `Species Category`)) +
geom_col(color = "white", size = 0.2) +
labs(x = "", y = "Proportion of Total Annual Biomass") +
coord_flip() +
scale_fill_gmri() +
scale_y_continuous(position = "right", ) +
guides("fill" = guide_legend(title.position = "top", title.hjust = 0.5)) +
theme(legend.position = "bottom")
Discrete Bodymass Bins
# Plot the cumulative plots by bodymass bins, by decade
decade_plots_size <- size_cumsums %>%
mutate(decade = floor_decade(Year)) %>%
split(.$decade) %>%
map(function(decade_group){
# Drop levels
levels_dropped <- decade_group %>%
mutate(Year = forcats:::fct_drop(Year))
totals_dropped <- filter(totals, Year %in% levels_dropped$Year) %>%
mutate(Year = forcats:::fct_drop(Year))
# Cumulative Biomass
biomass_plot <- levels_dropped %>%
ggplot(aes(x = `Weight Bin`, y = `Total Bin Biomass` / 1000)) +
geom_col(aes(fill = `Bodymass Bin (g)`)) +
geom_hline(data = totals_dropped, aes(yintercept = (`Total Annual Biomass` / 1000) / 2),
linetype = 2, alpha = 0.5, color = "royalblue") +
facet_wrap(~Year, ncol = 5) +
scale_y_continuous(labels = scales::comma_format()) +
scale_fill_gmri() +
theme(legend.position = "none") +
labs(y = "Total Bin Biomass (kg)")
# Proportion of Total Biomass
proportion_plot <- levels_dropped %>%
ggplot(aes(x = `Weight Bin`, y = `Percent of Annual Biomass`)) +
geom_col(aes(fill = `Bodymass Bin (g)`)) +
geom_hline(yintercept = 50, linetype = 2, alpha = 0.5, color = "royalblue") +
facet_wrap(~Year, ncol = 5) +
scale_fill_gmri() +
guides("color" = guide_legend(title.position = "top", title.hjust = 0.5)) +
theme(legend.position = "top")
stacked_plot <- biomass_plot / proportion_plot
return(list("mass" = biomass_plot,
"proportion" = proportion_plot))
})1980’s
# And plot
#decade_plots_size$`1980`$mass
decade_plots_size$`1980`$proportion
1990’s
#decade_plots_size$`1990`$mass
decade_plots_size$`1990`$proportion
2000’s
#decade_plots_size$`2000`$mass
decade_plots_size$`2000`$proportion
2010’s
#decade_plots_size$`2010`$mass
decade_plots_size$`2010`$proportion
Species Class Bins
#### Column Plots for the different classes
# And do the same for the species classes
decade_plots_class <- class_cumsums %>%
mutate(decade = floor_decade(Year)) %>%
split(.$decade) %>%
map(function(decade_group){
# Drop levels
levels_dropped <- decade_group %>%
mutate(Year = forcats:::fct_drop(Year))
totals_dropped <- filter(totals, Year %in% levels_dropped$Year) %>%
mutate(Year = forcats:::fct_drop(Year))
# Cumulative Biomass
biomass_plot <- levels_dropped %>%
ggplot(aes(x = `Species Class`, y = `Total Class Biomass` / 1000)) +
geom_col(aes(fill = `Species Category`)) +
geom_hline(data = totals_dropped, aes(yintercept = (`Total Annual Biomass` / 1000) / 2),
linetype = 2, alpha = 0.5, color = "royalblue") +
facet_wrap(~Year) +
scale_fill_gmri() +
scale_y_continuous(labels = scales::comma_format()) +
theme(legend.position = "none") +
labs(y = "Total Class Biomass (kg)")
# Proportion of Total Biomass
proportion_plot <- levels_dropped %>%
ggplot(aes(x = `Species Class`, y = `Percent of Annual Biomass`)) +
geom_col(aes(fill = `Species Category`)) +
geom_hline(yintercept = 50, linetype = 2, alpha = 0.5, color = "royalblue") +
facet_wrap(~Year) +
scale_fill_gmri() +
guides("color" = guide_legend(title.position = "top", title.hjust = 0.5)) +
theme(legend.position = "top")
stacked_plot <- biomass_plot / proportion_plot
return(list("mass" = biomass_plot,
"proportion" = proportion_plot))
})1980’s
# And plot
#decade_plots_class$`1980`$mass
decade_plots_class$`1980`$proportion
1990’s
#decade_plots_class$`1990`$mass
decade_plots_class$`1990`$proportion
2000’s
#decade_plots_class$`2000`$mass
decade_plots_class$`2000`$proportion
2010’s
#decade_plots_class$`2010`$mass
decade_plots_class$`2010`$proportion
Fish Length
#### Plotting CDP of lmin instead of mass bins ####
decade_plots_length <- lmin_cumsums %>%
mutate(decade = floor_decade(Year)) %>%
split(.$decade) %>%
map(function(decade_group){
# Drop levels
levels_dropped <- decade_group %>%
mutate(Year = forcats:::fct_drop(Year))
totals_dropped <- filter(totals, Year %in% levels_dropped$Year) %>%
mutate(Year = forcats:::fct_drop(Year))
# Cumulative Biomass
biomass_plot <- levels_dropped %>%
ggplot(aes(x = `Length (cm)`, y = `Cumulative Biomass (kg)`)) +
geom_step(group = 1, size = 1) +
geom_hline(data = totals_dropped, aes(yintercept = (`Total Annual Biomass` / 1000) / 2),
linetype = 2, alpha = 0.5, color = "royalblue") +
scale_y_continuous(labels = scales::comma_format()) +
facet_wrap(~Year, ncol = 5) +
theme(legend.position = "none")
# Proportion of Total Biomass
proportion_plot <- levels_dropped %>%
ggplot(aes(x = `Length (cm)`, y = `Proportion of Annual Biomass`)) +
geom_step(group = 1, size = 1) +
geom_hline(yintercept = 0.5, linetype = 2, alpha = 0.5, color = "royalblue") +
facet_wrap(~Year, ncol = 5)
stacked_plot <- biomass_plot / proportion_plot
return(list("mass" = biomass_plot,
"proportion" = proportion_plot))
})1980’s
# And plot
#decade_plots_length$`1980`$mass
decade_plots_length$`1980`$proportion
1990’s
#decade_plots_length$`1990`$mass
decade_plots_length$`1990`$proportion
2000’s
#decade_plots_length$`2000`$mass
decade_plots_length$`2000`$proportion
2010’s
#decade_plots_length$`2010`$mass
decade_plots_length$`2010`$proportion
Regional Temperature Relationship
For these we want to create polygons using the stratum for each region, that way there is an identical area used for calculation of regional anomalies.
OISST Timeseries for each region are then calculated in python, and available to access via box
# Format table to pair with temperature
results_5g <- min5g %>%
filter(Year != "all") %>%
mutate(
Year = as.numeric(Year),
region = case_when(
area == "all" ~ "All Regions",
area == "GB" ~ "Georges Bank",
area == "GoM" ~ "Gulf of Maine",
area == "MAB" ~ "Mid-Atlantic Bight",
area == "SNE" ~ "Southern New England"))
# Stratified Abundance Results
strat_5g <- min5g_strat %>%
filter(Year != "all") %>%
mutate(
Year = as.numeric(Year),
region = case_when(
area == "all" ~ "All Regions",
area == "GB" ~ "Georges Bank",
area == "GoM" ~ "Gulf of Maine",
area == "MAB" ~ "Mid-Atlantic Bight",
area == "SNE" ~ "Southern New England"))Temperature Timelines
# Format the temperature timeseries as well
regional_temps <- region_temps %>%
mutate(Year = lubridate::year(time),
est_month = lubridate::month(time),
season = case_when(
est_month %in% c(3:5) ~ "SPRING",
est_month %in% c(6:8) ~ "SUMMER",
est_month %in% c(9:11) ~ "FALL",
est_month %in% c(1,2,12) ~ "Winter" ))
# Plot regional timelines
region_temps %>%
ggplot(aes(time, sst)) +
geom_line(color = "gray50") +
facet_wrap(~region, ncol = 2) +
labs(x = "", y = "Sea Surface Temperature")
# Do some averaging to get spring/fall and annual temp and anoms
# and also some lagging action
# Spring primarily in March, April or May,
# Fall primarily in Sept, Oct, NovRegional SS Slopes & Temperature
# Format temperature data to combine with size spectrum groups
# Get the overall temperatures for the entire area
overall_temps <- regional_temps %>%
mutate(region = "All Regions") %>%
group_by(Year, region, season) %>%
summarise(avg_sst = mean(sst, na.rm = T),
avg_sst_clim = mean(sst_clim, na.rm = T),
avg_sst_anom = mean(sst_anom, na.rm = T)) %>%
ungroup()
# Get the regional and seasonal means
seasonal_temps <- regional_temps %>%
group_by(Year, region, season) %>%
summarise(
avg_sst = mean(sst, na.rm = T),
avg_sst_clim = mean(sst_clim, na.rm = T),
avg_sst_anom = mean(sst_anom, na.rm = T)) %>%
ungroup()
# Put together
seasonal_temps <- bind_rows(overall_temps, seasonal_temps) %>%
arrange(Year, region, season) %>%
mutate(season = str_to_title(season))
# Join temp to SS results
results_w_temps <- left_join(results_5g, seasonal_temps, by = c("Year", "region", "season")) %>%
drop_na(avg_sst)
# Join temp to stratified abundance results
strat_w_temps <- left_join(strat_5g, seasonal_temps, by = c("Year", "region", "season")) %>%
drop_na(avg_sst)Survey Abundances
These figures show how the size spectrum slopes vary using the observed abundances (numlen_adj) from the survey.
# Plot Survey Abundances
results_w_temps %>%
mutate(season = str_to_lower(season),
season = factor(season, levels = c("spring", "fall"))) %>%
ggplot(aes(avg_sst_anom, b)) +
geom_point() +
geom_smooth(formula = y ~ x,
size = 1,
method = "lm") +
stat_poly_eq(formula = y ~ x, parse = TRUE,
label.y = "bottom", label.x = "right") +
facet_grid(season ~ region) +
labs(x = "Seasonal Temperature Anomaly",
y = "Size Spectrum Slope (b)",
caption = "Size Spectrum Slopes Estimated from Observed Survey Catches")
Stratified Survey Abundances
# Plot Survey Abundances
strat_w_temps %>%
mutate(season = str_to_lower(season),
season = factor(season, levels = c("spring", "fall"))) %>%
ggplot(aes(avg_sst_anom, b)) +
geom_point() +
geom_smooth(formula = y ~ x,
size = 1,
method = "lm") +
stat_poly_eq(formula = y ~ x, parse = TRUE,
label.y = "bottom", label.x = "right") +
facet_grid(season ~ region) +
labs(x = "Seasonal Temperature Anomaly",
y = "Size Spectrum Slope (b)",
caption = "Size Spectrum Slopes Estimated from Stratified Abundances")
Regional SS Slope Timelines
These are the results from using a 5 gram bodysize cutoff.
Survey Abundances
These figures show how the size spectrum slopes vary using the observed abundances (numlen_adj) from the survey.
# # Plotting all 5gram cutoff results
min5g <- min5g %>%
mutate(Year = as.numeric(as.character(Year)),
season = case_when(season == "all" ~ "Spring + Fall",
season == "Spring" ~ "Spring",
season == "Fall" ~ "Fall"),
season = factor(season, levels = c("Spring", "Fall", "Spring + Fall")))
# # how the hline intercepts are created by facet grid
# filter(min5g, is.na(Year)) %>%
# split(.$season) %>%
# map(~ split(x = .x, f = .x$area))
min5g %>%
ggplot(aes(Year, b, color = area, shape = season)) +
geom_pointrange(aes(x = Year, y = b, ymin = confMin, ymax = confMax), alpha = 0.5) +
geom_hline(data = filter(min5g, is.na(Year)),
aes(yintercept = b), linetype = 2, alpha = 0.7, size = 1) +
facet_grid(season ~ area) +
geom_smooth(method = "gam", show.legend = FALSE, se = FALSE,
formula = y ~ s(x, bs = "cs")) +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, size= 6, vjust = 0.5)) +
scale_color_gmri() +
labs() +
guides(shape = "none",
color = "none") +
labs(x = "",
caption = "Horizontal line represents the mean across all years for facet group.",
y = "Size Spectrum Slope (b)") 
Stratified Survey Abundances
These figures show how the size spectrum slopes vary using the area-stratified abundances. These scale out catch annual catch rates from the survey strata to the total areas.
# # Plotting all 5gram cutoff results
min5g_strat <- min5g_strat %>%
mutate(Year = as.numeric(as.character(Year)),
season = case_when(season == "all" ~ "Spring + Fall",
season == "Spring" ~ "Spring",
season == "Fall" ~ "Fall"),
season = factor(season, levels = c("Spring", "Fall", "Spring + Fall")))
min5g_strat %>%
ggplot(aes(Year, b, color = area, shape = season)) +
geom_pointrange(aes(x = Year, y = b, ymin = confMin, ymax = confMax), alpha = 0.5) +
geom_hline(data = filter(min5g, is.na(Year)),
aes(yintercept = b), linetype = 2, alpha = 0.7, size = 1) +
facet_grid(season ~ area) +
geom_smooth(method = "gam", show.legend = FALSE, se = FALSE,
formula = y ~ s(x, bs = "cs")) +
theme_bw() +
theme(axis.text.x = element_text(angle = 90, size= 6, vjust = 0.5)) +
scale_color_gmri() +
labs() +
guides(shape = "none",
color = "none") +
labs(x = "",
caption = "Horizontal line represents the mean across all years for facet group.",
y = "Size Spectrum Slope (b)") 